Esta es una guía detallada de como poder realizar el trabajo práctico N°1 del curso Planificación Territorial.
El objetivo es … jardines … ciudad en 15 minutos … ver cuanto se demora a pie.
1 Organización
Importante
Este informe esta pensado en ejecutarse dentro de una carpeta que contenga el trabajo (por ejemplo, trabajo1) y DENTRO DE ESTA una carpeta llamada datos, la cual contendrá todas las capas para el análisis.
A modo de ejemplo, el trabajo debería tener un orden así
Paquete para el manejo de datos. Personalmente prefiero cargar tidyverse ya que incluye otras librerías importantes como lubridate (para manejo de fechas) o ggplot2, para visualización.
2
Librería para analizar teoría de grafos. Encontré una guía introductoria (pero en inglés) aquí.
3
Librería que maneja los datos espaciales. Piensen los datos vectoriales como una tabla excel con atributos (como Nombre, Descripción, Region, etc.) y que al final tiene una columna de geometría (usualmente llamada .geo o geometry), la que provee a un dato de una “forma”.
4
Paquete para hacer los hexágonos.
5
Librería para manejar datos ráster.
6
Para obtener datos de elevación y de Open Street Map (OSM).
7
Esta es para hacer cartografía, tal como se hace en QGIS. Pero con esta librería se hacen los mapas interactivos también.
Nota
Se puede hacer click en los números para una “selección” visual.
2.2.2 Situar el Ambiente
Antes de comenzar, dentro del script deberemos setear el ambiente hacia donde están los datos. Revisar Sección 1 para seguir el siguiente ejemplo.
Si el script (scprit.r) se encuentra dentro de la carpeta trabajo1, pero no dentro de la carpeta datos
setwd("datos/")
Si el script no se encuentra dentro de la carpeta trabajo1, se deberá hacer alusión a la ruta completa de los datos. En mi caso:
setwd("C:/uni/oit/t1/datos/")
Si el script se encuentra dentro de la carpeta datos, en teoría no debería hacerse nada.
Advertencia
Windows indica las rutas de las carpetas con el backslash (\) pero R nos los maneja bien. Para indicar la ruta de una carpeta se debe cambiar este símbolo por un slash (/).
2.2.3 Cargar las capas
Una vez seteados en la carpeta correcta, podemos cargar todas las capas y trabajar en base a eso. Para los datos .shp se deberá ocupar st_read() desde el paquete sf y para datos raster la función rast() de terra.
En el caso de las capas de LUC y Unidad Vecinal de Valdivia, que se encuentran en formato .gpkg, se deberá indicar un argumento extra dentro de st_read(). Esto ya que el formato geopackage es un contenedor que puede almacenar múltiples capas y otros datos relacionados, todo en un único archivo. Es por esto que hay que indicar que capa se quiere acceder. Para ver el listado de capas que tiene un .gpkg, se puede usar la funcion st_layers():
st_layers("LUC_Valdivia.gpkg")
Driver: GPKG
Available layers:
layer_name geometry_type features fields crs_name
1 luc_valdivia Multi Polygon 1 12 WGS 84
Podemos ver que el archivo LUC_Valdivia.gpkg tiene solo una capa, llamada luc_valdivia (a través del campo layer_name). En el caso de las unidades vecinales
st_layers("Unidad_vecinal_Valdi.gpkg")
Driver: GPKG
Available layers:
layer_name geometry_type features fields crs_name
1 unidad_vecinal 3D Multi Polygon 57 11 WGS 84
la capa se llama unidad_vecinal. Con esta información podemos cargar correctamente los datos:
Para poder hacer el informe, debemos tener en consideración la proyección de las capas. Para que el flujo funcione, TODAS las capas deben estar en EPSG:4326, el cual corresponder a la proyección mundial WGS 84, el cual se compone de latitud y longitud (en grados).
Datos Vectoriales
Para ver el CRS de un dato vectorial se podría ocupar la función st_crs(), el cual mostrará toda la información de la proyección.
st_crs(uvVald)
Coordinate Reference System:
User input: WGS 84
wkt:
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
MEMBER["World Geodetic System 1984 (G2296)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]
Sin embargo, como esto muestra toda la información, accederemos únicamente al EPSG accediendo con el signo peso ($).
st_crs(uvVald)$epsg
[1] 4326
Esta capa ya se encuentra en el sistema de coordenadas. De querer confirmar esto, se puede comparar con el operador de igualdad (==). Este operador compara si dos valores son iguales. Devuelve TRUE si son iguales y FALSE si son distintos.
st_crs(uvVald)$epsg ==4326
[1] TRUE
Esto confirma que la capa se encuentra en el sistema de proyección correspondiente. Veamos las otras capas
st_crs(lucVald)$epsg ==4326
[1] TRUE
st_crs(jardinJunji)$epsg ==4326
[1] TRUE
st_crs(jardinIntegra)$epsg ==4326
[1] TRUE
Esto indica que todas las capas se encuentran en WGS 84.
Dato Ráster
De la misma manera, terra tiene la función crs() que con el argumento describe = TRUE y accediendo al código (code) con $ podemos comparar de la misma manera
crs(dem, describe = T)$code ==4326
[1] TRUE
Con esto confirmamos que todas las capas están en la misma proyección.
2.2.5 Extensión
El objetivo del trabajo está pensado en las unidades vecinales de Valdivia (uvVald) pero dentro de los límites de la ciudad (lucVald). Visualicemos como se situan todas las capas.
Make this Notebook Trusted to load map: File -> Trust Notebook
A primera vista podemos dos desafíos:
La enorme cantidad de establecimientos en jardines infantiles. Esto ocurre ya que ambas campas (Junji e Integra) se encuentran a nivel nacional.
También, el LUC nos provee de los límites a trabajar mientras que las unidades vecinales (que son más grandes), proporciona las “particiones” que tendría la ciudad de Valdivia.
Para ambos problemas se pueden ocupar operaciones espaciales, como por ejemplo, el corte o la superposición. El flujo sería:
Cortar las Unidades Vecinales con las delimitaciones de Valdivia (dadas por el LUC).
A través del LUC de Valdivia, cortar ambas capas de jardines infantiles (Junji e Integra).
De este modo, tendríamos todas las capas con las extensiones correspondientes.